The SIR model in pomp

Licensed under the Creative Commons attribution-noncommercial license, http://creativecommons.org/licenses/by-nc/3.0/. Please share and remix noncommercially, mentioning its origin.
Produced using pomp version 0.71.1.

Compartmental models in epidemiology

The SIR model

Below is a diagram of the SIR model. The host population is divided into three classes according to their infection status: S, susceptible hosts; I, infected (and infectious) hosts; R, recovered and immune hosts. The S\(\to\)I rate, \(\lambda\), called the force of infection, depends on the number of infectious individuals according to \(\lambda(t)={\beta\,I}/{N}\). The I\(\to\)R, or recovery, rate is \(\gamma\). The case reports, \(C\), result from a process by which infections are recorded with probability \(\rho\). Since diagnosed cases are treated with bed-rest and hence removed, infections are counted upon transition to R.

Implementing compartmental models in pomp

Load the data

baseurl <- "http://kinglab.eeb.lsa.umich.edu/SBIED/"
url <- paste0(baseurl,"data/bs_flu.csv")
read.csv(url,comment.char="#",colClasses=c(date="Date")) %>%
  subset(status=="confined") %>%
  mutate(day=as.integer(date-min(date))) %>%
  subset(select=-c(status,date)) -> bbs
bbs %>% ggplot(aes(x=day,y=cases))+geom_line()+geom_point()

Set up the process model. We need a model that will simulate the process from time \(t\) to time \(t+{{\Delta}{t}}\).

sir_step <- "
  double t1 = rbinom(S,1-exp(-beta*I/pop*dt));
  double t2 = rbinom(I,1-exp(-gamma*dt));
  S -= t1;
  I += t1 - t2;
  R += t2;
"

Let’s assume that the data represent incidence, i.e., the number of new infections occurring on a given day.

sir_step <- "
  double t1 = rbinom(S,1-exp(-Beta*I/763*dt));
  double t2 = rbinom(I,1-exp(-gamma*dt));
  S -= t1;
  I += t1 - t2;
  R += t2;
  H += t1;
"

To initialize the state process, we write

sir_init <- "
  S = 762;
  I = 1;
  R = 0;
  H = 0;
"

We add these to the data to make a pomp object:

bbs %<>% pomp(times="day",t0=0,
              rprocess=euler.sim(step.fun=Csnippet(sir_step),delta.t=0.1),
              initializer=Csnippet(sir_init), zeronames="H",
              statenames=c("S","I","R","H"),
              paramnames=c("Beta","gamma"))

H now accumulates the new infections. The incidence on day \(t\) is \(H(t+1)-H(t)\).

We’ll model the data with some degree, \(\rho\) of under-reporting: \[\text{cases}_t \sim {\mathrm{Binomial}\left(H(t+1)-H(t),\rho\right)}.\]

As before, we must write both a dmeasure and an rmeasure component:

dmeas <- "lik = dbinom(cases,H,rho,give_log);"
rmeas <- "cases = rbinom(H,rho);"

and put these into bbs:

bbs %<>% pomp(rmeasure=Csnippet(rmeas),dmeasure=Csnippet(dmeas),
              statenames="H",paramnames="rho")

simulate(bbs,params=c(Beta=2,gamma=0.1,rho=0.9),nsim=20,as=TRUE,include=TRUE) %>%
  ggplot(aes(x=time,y=cases,group=sim,color=sim=="data"))+
  geom_line()+guides(color=FALSE)

The parameters are constrained to be positive, and \(\rho < 1\). We’ll find it useful to transform the parameters onto a scale on which there are no such constraints. The following accomplish this.

toEst <- "
 TBeta = log(Beta);
 Tgamma = log(gamma);
 Trho = logit(rho);
"

fromEst <- "
 TBeta = exp(Beta);
 Tgamma = exp(gamma);
 Trho = expit(rho);
"

bbs %<>% pomp(toEstimationScale=Csnippet(toEst),
              fromEstimationScale=Csnippet(fromEst),
              paramnames=c("Beta","gamma","rho"))

Sequential Monte Carlo in pomp

In pomp, the basic particle filter is implemented in the command pfilter. We must choose the number of particles to use by setting the Np argument.

pf <- pfilter(bbs,Np=1000,params=c(Beta=2,gamma=0.1,rho=0.9))
logLik(pf)

## [1] -381.0412

We can run a few particle filters to get an estimate of the Monte Carlo variability:

pf <- replicate(n=10,pfilter(bbs,Np=1000,params=c(Beta=2,gamma=0.1,rho=0.9)))
ll <- sapply(pf,logLik)
logmeanexp(ll,se=TRUE)

##                      se 
## -380.709774    1.257268

Note that we’re careful here to counteract Jensen’s inequality. The particle filter gives us an unbiased estimate of the likelihood, not of the log-likelihood.

To get an idea of what the likelihood surface looks like in the neighborhood of the default parameter set supplied by bbs, we can construct some likelihood slices. We’ll make slices in the \(\beta\) and \(\gamma\) directions. Both slices will pass through the default parameter set.

require(pomp)
require(foreach)
require(doMC)

sliceDesign(
  c(Beta=2,gamma=0.1,rho=0.9),
  Beta=rep(seq(from=0.5,to=1,length=40),each=3),
  gamma=rep(seq(from=0.01,to=0.25,length=40),each=3)) -> p

registerDoMC(cores=4)        ## number of cores on your machine

set.seed(998468235L,kind="L'Ecuyer")
mcopts <- list(preschedule=FALSE,set.seed=TRUE)

foreach (theta=iter(p,"row"),.combine=rbind,
         .inorder=FALSE,.options.multicore=mcopts) %dopar% 
  {
    pfilter(bbs,params=unlist(theta),Np=1000) -> pf
    theta$loglik <- logLik(pf)
    theta
    } -> p

p %>%
  mutate(slice=as.character(slice)) %>%
  melt(id=c("slice","loglik")) %>% 
  subset(variable==slice) %>% 
  dlply(~slice, function (df) {
    ggplot(data=df,mapping=aes(x=value,y=loglik))+
      geom_point()+
      geom_smooth()+
      labs(x=unique(df$variable))
  })

Exercise

Add likelihood slices along the \(\rho\) direction.

Clearly, the default parameter set is not particularly close to the MLE. One way to find the MLE is to try optimizing the estimated likelihood directly. There are many optimization algorithms to choose from, and many implemented in R.

Two issues arise immediately. First, the particle filter gives us a stochastic estimate of the likelihood. We can reduce this variability by making Np larger, but we cannot make it go away. If we use a deterministic optimizer (i.e., one that assumes the objective function is evaluated deterministically), then we must control this variability somehow. For example, we can fix the seed of the pseudo-random number generator. A side effect will be that the objective function becomes jagged, marked by many small local knolls and pits. Alternatively, we can use a stochastic optimization algorithm, with which we will be only be able to obtain estimates of our MLE. This is the trade-off between a noisy and a rough objective function. Second, because the particle filter gives us just an estimate of the likelihood and no information about the derivative, we must choose an algorithm that is “derivative-free”. There are many such, but we can expect less efficiency than would be possible with derivative information. Note that finite differencing is not an especially promising way of constructing derivatives. The price would be a \(n\)-fold increase in cpu time, where \(n\) is the dimension of the parameter space. Also, since the likelihood is only estimated, we would expect the derivative estimates to be noisy.

Here, let’s opt for deterministic optimization of a rough function. We’ll try using optim’s default method: Nelder-Mead.

coef(bbs) <- c(Beta=2,gamma=0.1,rho=0.9)

neg.ll <- function (par, est, ...) {
  ## parameters to be estimated are named in 'est'
  allpars <- coef(bbs,transform=TRUE)
  allpars[est] <- par
  try(
    pfilter(
      bbs,
      params=partrans(bbs,allpars,dir="fromEst"),
      ...
      )
    ) -> pf
  if (inherits(pf,"try-error")) {
    1e10 ## a big number
    } else {
      -logLik(pf)
      }
  }

require(plyr)
## use Nelder-Mead with fixed RNG seed
fit <- optim(
  par=c(-1.1, 0.33, 2.2),
  est=c("gamma","Beta","rho"),
  Np=200,
  fn=neg.ll,
  seed=3488755L,
  method="Nelder-Mead",
  control=list(maxit=400,trace=0)
  )

mle <- bbs
coef(mle,c("gamma","Beta","rho"),trans=T) <- fit$par
coef(mle,c("gamma","Beta","rho")) ## point estimate

##     gamma      Beta       rho 
## 0.4016605 1.3944346 0.9038721

mle %>% 
  simulate(nsim=9,as.data.frame=TRUE,include=TRUE) -> sims

lls <- replicate(n=5,logLik(pfilter(mle,Np=1000)))
logmeanexp(lls,se=TRUE)

##                        se 
## -372.2350171    0.3285705

Some simulations at these parameters are shown next:

sims %>%
  ggplot(aes(x=time,y=cases,group=sim,color=sim=="data"))+
  geom_line()+
  geom_text(x=9,y=320,
            label=paste("log~L==",round(ll[1],1),"%+-%",round(ll[2],1)),
            color="black",
            parse=T,hjust=0)

Exercise

Use simulated annealing to maximize the likelihood. Be sure to try several starting guesses.

Exercise

Construct likelihood slices on \(\beta\) and \(\gamma\) through the MLE you found above.

References